On the Saturation of the Nonlinear Refractive Index in Atomic Gases 



C. Kohler, 1 Roland Guichard, 2 Emmanuel Lorin, 3 Szczepan Chelkowski, 4 
Andre D. Bandrauk, Luc Berge, 1 and Stefan Skupin 5 ' 6 

1 CEA-DAM, DIF, F-91297 Arpajon, France 
2 CNRS, UMR 7614, LCPMR, 75231 Pans Cedex 05, France 
Carleton University, Ottawa (Ontario) K1S 5B6, Canada 
' Departement de chimie, Universite de Sherbrooke, Sherbrooke (Quebec) J1K 2R1, Canada 
''Max Planck Institute for the Physics of Complex Systems, 01187 Dresden, Germany 
6 Friedrich Schiller University, Institute of Condensed Matter Theory and Optics, 07743 Jena, Germany 

Motivated by the ongoing controversy on the origin of the nonlinear index saturation and sub- 
sequent intensity clamping in femtosecond filaments, we study the atomic nonlinear polarization 
induced by a high-intensity and ultrashort laser pulse in hydrogen by numerically solving the time 
dependent Schrodinger equation. Special emphasis is given to the efficient modeling of the nonlinear 
polarization at central laser frequency corresponding to 800 nm wavelength. Here, the recently pro- 
posed model of the Higher-Order Kerr Effect (HOKE) and two versions of the Standard model for 
femtosecond filamentation, including either a multi-photon or tunnel ionization rate, are compared. 
We find that around the clamping intensity the instantaneous HOKE model does not reproduce 
the temporal structure of the nonlinear response obtained from the quantum mechanical results. 
In contrast, the non-instantaneous charge contributions included in the Standard models ensure a 
reasonable quantitative agreement. Therefore, the physical origin for the observed saturation of 
the overall electron response is confirmed to mainly result from contributions of free or nearly free 
electrons. 



PACS numbers: 32.80.-t,51.50.+v 

I. INTRODUCTION 

Accurate modeling of high intensity laser pulse prop- 
agation in a medium requires the accurate knowledge of 
its optical response. For laser intensities well below the 
atomic ionization threshold this response is given by the 
induced polarization involving bound electrons, usually 
modeled by a Taylor expansion in terms of the driv- 
ing electric field. Restriction of the expansion to the 
third order nonlinear term for isotropic media leads to 
the well known intensity-dependent Kerr refractive in- 
dex change PQ. This simplest approach for the polariza- 
tion already provides all key ingredients to describe fun- 
damental nonlinear phenomena like modulational insta- 
bility, self- focusing and self-phase modulation 2J, which 
manifest in soliton formation, pulse compression or self- 
similar collapse. The latter effect, however, demonstrates 
the inherent incompleteness of this simplest approach as 
it predicts an unphysical blow-up of the beam intensity, 
which is intimately linked to the onset of laser filamen- 
tation. 

A theoretical description of laser filamentation there- 
fore requires an accurate modeling of the saturation of 
the optical nonlinearity, i.e., the induced variation of 
the nonlinear refractive index with respect to the laser 
intensity. Originally achieved by accounting for photo- 
ionization and subsequent interaction of the laser pulse 
with a free electron plasma [3HSj , this so-called Standard 
model for nonlinear refractive index saturation was re- 
cently challenged by the measurement of higher order 



Kerr coefficients (TUHH]. Being negatively valued, these 
coefficients could give rise to an alternative explanation 
for the saturation of the optical nonlinearity. This strik- 
ingly simple approach, based on higher order contribu- 
tions in the Taylor expansion of the polarization in terms 
of the driving electric field, is called the Higher-Order 
Kerr Effect (HOKE) model. Besides the apparent dif- 
ferences in the model equations, there is a fundamental 
physical difference between HOKE and Standard model. 
While in the latter model a free electron plasma is respon- 
sible for the arrest of the self-focusing, HOKE terms orig- 
inate exclusively from the response of bound electrons. 
Thus, the HOKE model questions the very nature of 
the previous understanding of femtosecond filamentation 
as a dynamical balance between optical Kerr effect and 
self-generated electron plasma [T2TI20) . From the exper- 
imental side, numerous attempts have been undertaken 
to decide in favor for either bound or ionized electron 
contributions 2f -28j without common agreement. Re- 
cent progress was achieved by direct numerical solution of 
the time dependent Schrodinger equation (TDSE), which 
showed the excitation of atoms into highly excited Ryd- 
berg and Kramers-Henneberger states. Although bound, 
these states lie at the frontier above which excited elec- 
trons may be fully freed, and they behave very similarly 
to states in the continuum [T31 Q23 HO] • 

Here, we propose a direct confrontation of both HOKE 
and Standard models with rigorous quantum mechanical 
calculations. To this end, we investigate the saturation 
of the nonlinear polarization by numerically solving the 
TDSE for an atomic hydrogen model, for the sake of sim- 
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plicity first in one dimension (ID), and later in full three 
dimensions (3D) for confirmation. It turns out that a 
consistent definition of the HOKE coefficients is difficult, 
because the validity of the Taylor expansion of the non- 
linear response in terms of the electric field breaks down 
at intensities well below the clamping one. Moreover, 
the instantaneous character of the HOKE model contra- 
dicts the quantum mechanical results. In contrast, we 
achieve reasonable agreement when applying the Stan- 
dard model with multi-photon or tunnel ionization rate. 
Although the exact definition of bound and ionized elec- 
trons remains ambiguous for atoms in strong laser fields, 
the saturating contribution to the nonlinear refractive in- 
dex is still well described by a simple Drude plasma model 
assuming continuum electrons as encountered once the 
pulse has passed the interaction region. 

The paper is organized as follows. In Sec. [Hj we intro- 
duce our TDSE modeling of the nonlinear polarization as 
well as the phenomenological models currently discussed 
in the literature. For sake of simplicity, we focus on ID 
TDSE description of atomic hydrogen. Then, we con- 
front the phenomenological models (HOKE and Standard 
model) with rigorous TDSE results in Sec. Ill Finally, 



we confirm that our findings hold for more realistic 3D 
TDSE calculations in Sec. EI 



II. MODELING NONLINEAR POLARIZATION 

A. ID Quantum Mechanical Treatment 

Let us start with ID quantum mechanical calculations. 
We assume that the single non-relativistic electron of the 
ID hydrogen model in an external laser E(t) is described 
by its wave function ^(x,t), governed by the TDSE 



(1) 



Here Hq = —h 2 d 2 /(2m e ) + V(x) with electron mass m e , 
and the interaction with the external laser field in veloc- 
ity gauge is formulated in the dipole approximation, i.e. 
Hint = ihq e A(t)d x /m e involving the electron charge q e 
and the magnetic potential A(t) — — J_ E(t')dt' . We 

use a soft-core potential V(x) = — g 2 /(47re \/a: 2 + 2a 2 ,) 
adjusted to match the ionization energy of hydrogen, 
where ao is the Bohr radius. 

The macroscopic polarization is then given by P(t) = 
p at q e (^(t)\x\^!(t)), where p at is the atomic density at 
ambient pressure. 

As proposed in [IS], we distinguish between bound 
and continuum electrons by projecting the electronic 
wavefunction ^ on bound eigenfunctions of the field free 
Hamiltonian Ho\ip n ) = E n \ip n ), such that 

|*B>=5>n><W*>, |*c) = |*> - |*B>, (2) 

n 

where we used states up to principal quantum number 
n = 52. Then, the macroscopic polarization P(t) can 



be decomposed into bound, bound-continuum and con- 
tinuum part, P(t) = P BB {t) + 29?[P BC (t)] + Pcc(t), ac- 
cording to Eq. pj). For example, 

Pcc{t)=PaM*c{t)\x\* C {t))- (3) 

To observe the dependence of P(t) on the laser pulse 
intensity, we numerically integrate Eq. ([!]) for Gaussian 
pulses 



E(t) = E Q e 



-{t/rf 



(w i) 



(4) 



with central frequency ujq = 2t:c/Xq (Ao = 800 nm, i.e. 
cycle period of 2.67 fs), and peak intensities I ~ |Po| 2 
varying from 10~ 3 TW/cm 2 to 70 TW/cm 2 for pulse 
durations r = 10 fs, 20 fs and 100 fs. This intensity range 
ensures partial ionization by multi-photon and/or tunnel 
processes as complete ionization for hydrogen appears at 
/ = 120 TW/cm 2 by over-the-barrier ionization [2"5] . 

In order to access the macroscopic nonlinear polariza- 
tion Pnl, we first extract the linear response Phin{t) = 
CqX^ E(t), where we neglect linear dispersion and ab- 
sorption, from low intensity calculations I < 1 TW/cm 2 
yielding x^ 1 -* = 6.88 x 10~ 4 . Then, the nonlinear part of 
the polarization can be obtained with P/vx(i) = Pit) — 
eox^E(t). For the choice of the carrier envelope phase 
offset in our test pulses [Eq. Q], $1[Pnl (wo)] is responsi- 
ble for self-focusing (?R.[Pnl (wq)] > 0) or self-defocusing 
(5R[PjV£ ( w o)] < 0) action, i.e., a positive (resp. negative) 
nonlinear induced refractive index change. ^[Pnl (wo)] 
describes nonlinear absorption, which is weak and will be 
neglected throughout this paper. 



B. Phenomenological Models 

Let us now investigate the possibility to mimic the be- 
havior of the nonlinear polarization from TDSE calcula- 
tions by means of approximate models. Such phenomeno- 
logical description of the nonlinear response is extremely 
important for the efficient numerical modeling of high 
intense laser pulse propagation. As already mentioned 
above, we will approach the TDSE results in terms of 
two controversial models, currently used in the context 
of femtosecond filamentation. These models involve only 
a single set of few parameters for all pulse durations and 
intensities. First, the HOKE model 

Pnl® = e QX {3} E(tf + e oX (5) E(t) 5 + e oX (7) E(t) 7 + ■■■ 

= e J2 X U) E(ty (5) 

3=3,5,- 

is fully determined by the nonlinear coefficients ■ Sec- 
ond, the Standard model 

Pnl® = e oX (3) E(t) 3 - q 2 J (m eW 2 ) p(t)E(t), (6) 

requires an additional equation to describe the free elec- 
tron density 



d tP (t) = w(t) (pat ~ P(t)) ■ 



(7) 
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Here, we will either apply a multi-photon ionization 
(MPI) rate [30] 



= WMPi(I(t)) = a K I(t) 



K 



(8) 



with photon number K and ionization cross section <jk 
or a field-dependent tunnel ionization rate [311 132j 

w = w tunnel {E(t)) = a/\E{t)\ exp {-p/\E{t)\) (9) 

involving two parameters a, j3. 



Extraction of Phenomenological Model 
Parameters from ID TDSE Results 



In order to confront the phenomenological models with 
TDSE results, we first extract the susceptibilities x^ f° r 
the HOKE model and additionally the MPI cross section 
ok resp. the tunnel ionization parameters a, f3 for the 
two Standard models from the quantum mechanical data. 

For the extraction of the x we use the assumption 
that different orders P^|(i) = e X (j) E(t) j of the nonlin- 
ear polarization become important subsequently for in- 
creasing field strength. Additionally, we assume that the 
characteristic peak at the jth harmonic ju>o in P^ L is 
mainly caused by x ■ Then, one can deduce consecu- 
tively the x using 



(10) 



valid for the Gaussian pulses of Eq. Q. Alternatively, 
we can evaluate Pjvl(wo, I) at different peak intensities I 
and infer the values of %W by solving a system of linear 
equations 

P NL (uo, h) = e FT{ X i3) Ef(t) + x^E^t) + • • ■ }U, 
Pnl(uo, h) = e FT{ X {3) El(t) + X (5) E 5 2 (t) + • • • }\ Uo , 
••. = ..., (11) 

where FT{} denotes the Fourier transform. Note, that 
the former approach evaluates the susceptibilities at dif- 
ferent harmonic frequencies jwo, whereas the latter one 
gives access to all x directly at u)q. For purely instan- 
taneous HOKE terms of Eq. ([5| , both approaches should 
yield the same X ■ 

Let us first evaluate x a t the harmonic frequencies 



jujq. Results are presented in Figs. E](a)-(c), where x 



(3) 



X^ and x^ 7 ' are plotted for different pulse durations r 
vs. peak intensity. We can observe, how difficult it is 
to extract meaningful coefficients x being independent 
of intensity and pulse duration in the range where Pnl 
becomes relevant, i.e. I > 1 TW/cm 2 . The fact that 
it is not possible to describe the nonlinear polarization 
in terms of a Taylor expansion at intensities above 10 
TW/cm 2 already indicates significant changes in the elec- 
tronic configuration of the atomic system. This involves, 



e.g., the enhanced population of Stark shifted Rydberg 
states, staying close to the ionization threshold due to a 
comparable increase of the ponderomotive potential [53] . 
Moreover, the observed dependency of the higher order 
X^' on the pulse duration reflects strong nonlinear dis- 
persion, i.e., a noninstantaneous character of the nonlin- 
ear response at higher intensities. Due to this nonlinear 
dispersion, and because we are searching for a descrip- 
tion of Pnl at ujq, we resort to the latter approach of 
Eqs. (11) in the following. By doing so, we can only 



extract a minimum set of consistent parameters, namely 



X (3) = 2.3 x 10~ 25 (m/V) 2 , 
X (5) = -2.5 x 10~ 45 (m/V) 4 



(12) 
(13) 



all higher order x being discarded, since all of them 
strongly fluctuate along the intensity range. We will 
further refer to this simple model as HOKE_35. Note 
that we will later introduce another set of parameters, 
HOKE_FIT, which includes nonlinear coefficients up to 
X^ 9 \ chosen to obtain a more accurate description of the 
nonlinear saturation. 

For the Standard model with MPI rate, in addition to 
X^ 3 \ the photon number K and the cross section ax has 
to be determined from the TDSE simulation data. For 
pulses with wavelength Ao = 800 nm and the ionization 
potential of hydrogen of Uh = 13.6 eV one obtains a min- 
imum photon number of K = mod (Uul^Ja + 1) = 9. 
However, we find the least dependence of <jk on the peak 
intensity when employing a reduced photon number of 
K = 6. This is plausible, since for the here employed 
intensities the Keldysh parameter 7 is in the order of 
one, indicating a transient regime between MPI (7 3> 1) 
and tunnel ionization (7 <C 1), which leads to a reduced 
exponent in the intensity dependence of the ionization 
probability. In addition, present atomic resonances con- 
necting the ground and first excited state are of similar 
order (five and seven photon resonance in ID and 3D, 
respectively) (34] [35] . Because, we are only interested in 
the best (phenomenological) description of TDSE results, 
we will employ the MPI rate with K — 6 in the following. 
To calculate a^, we match the free electron density at the 



end of the pulse (t = 
with the TDSE one: 



tend) obtained with the MPI rate 



p(t -> t end ) = / dx\^ c (x,t end 



(14) 



where ^c(x,i) is the continuum part of the wavefunc- 
tion. Interestingly, the results for all three pulse dura- 
tions agree, up to strong modulations of <jq at increasing 
intensity [Fig. [11(d)]. We attribute these modulations to 
resonant excitation of high lying Rydberg states, where 
the energy of an integer number of photons matches the 
energy gap between the ground and a (Stark shifted) Ry- 
dberg state [20] . We thus select an average value for the 
Standard model, 



cr 6 = 8.6 x 10- di (cm 2 /TW)7s, 



(15) 
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FIG. 1: Model parameters extracted from ID TDSE simulations, (a) x {3 \ (b) X (5) , (c) X (7) for r = 10 fs (blue), r = 20 fs (red) 
and r = 100 fs (green) as function of the peak intensity /. (d) The MPI cross section ag extracted from ID TDSE simulations 
[same color coding as in (a)-(c)]. 



matching cr 6 at I — 29 TW/cm 2 . Our choice offers a 
reasonably good fit of the nonlinear polarization close to 
the clamping intensity, where Pnl changes sign. Be- 
low this intensity charge contributions are negligible, 
thus deviating erg do not matter. For higher intensi- 
ties we will overestimate the charge, however, at least 
in the context of femtosecond filamentation, this regime 
is expected to be barely relevant as well. In any case, 
complete ionization of the hydrogen atom is expected at 
I > 1.2 x 10 14 W/cm 2 by over-the-barrier ionization [2"9"| . 

To extract the tunnel ionization parameters we calcu- 
late the free electron density at the end of the pulse ac- 
cording to the tunnel ionization rate for all peak intensity 
values employed in the TDSE simulations and compare 
p(t t end ) to the TDSE data for r = 20 fs. Then, a 
and (3 are optimized such that the sum of the squared dif- 
ferences from each comparison is minimized. Performing 
this least squares fit, we extract 

a = 8.9 x 10 26 V/ms, (16) 
/3 = 1.6 x 10 n V/m, (17) 

for the parameters entering Eq. These values are 

used for all pulse durations in the following. 



III. CONFRONTATION OF 
PHENOMENOLOGICAL MODELS WITH ID 
TDSE RESULTS 

A. Behavior in Fourier Domain 

Having all necessary parameters at hand, let us now 
confront the HOKE and plasma-based (Standard) de- 
scription of the nonlinear polarization with TDSE results. 
In Figs. J2](a) and (b), we compare the nonlinear polar- 
ization P/vl(wo) at center frequency of the laser obtained 
from TDSE simulations with the approximate models for 
t = 20 fs and r = 100 fs, respectively. Shorter pulses 
(t = 10 fs) yield similar behavior (not shown). For the 
TDSE results with t = 20 fs (solid black line) Pnl(.ijJo) is 
positive and increases (focusing action) with peak inten- 
sity up to / ~ 25 TW/cm 2 . Then, saturation and final 
change of sign occur at I sat ~ 30 TW/cm 2 , signaling a 
qualitative change from focusing to defocusing behavior. 
This generic behavior also holds for r = 100 fs with de- 
creased [J sat ~ 20 TW/cm 2 , Fig. |] (b)] and r = 10 fs 
with increased [I sat ~ 40 TW/cm 2 , not shown] saturation 
intensity values. Moreover, we observe oscillatory struc- 
tures of P N l(ujo) at I ~ 35 TW/cm 2 (r = 20 fs) and 
I ~ 22 TW/cm 2 (t = 100 fs), which can be attributed 
to decreased ionization probabilities due to channel clo- 
sure [35]. The HOKE_35 model, defined by extracted 
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FIG. 2: ID simulation results. For (a) r = 20 fs and (b) r = 100 fs the nonlinear polarization at ljo is shown for TDSE 
calculations (solid black lines with circles) and compared to the Standard models (solid red lines = MPI rate, solid blue = 
tunnel rate), the HOKE_35 model (dotted green lines) containing a x' 3 ' and contribution only and the fitted HOKE_FIT 
model (dash-dotted green lines), (c) distinguishes bound (green), continuum (blue) and bound-continuum (red) electronic 
contributions to the nonlinear polarization (black) for r = 20 fs (see text for details). Additionally, the continuum contribution 
from the Standard model with tunnel ionization rate is shown (dashed blue line). In (d) the free electron density according 
to TDSE simulations (black) is compared to the one obtained from the MPI rate (red) and tunnel rate (blue) as used in the 
Standard models. 



X^ 3 ' and x^ 5 ' coefficients [see Eqs. (12), (13)] (green dot 



ted line) yields significantly lower intensity thresholds for 
saturation and sign inversion of Pnl(s^o) than observed 
in the TDSE calculations. Because this discrepancy ren- 
ders the HOKE_35 unsuitable for practical applications, 
we suggest a different (purely phenomenological) set of 
parameters HOKEJTT 



X (3) = 2.3 x 1(T 25 (m/V) 2 , 
X {5) = -2.0 x lCT 45 (m/V) 4 
X (7) = 5.3 x 1CT 65 (m/V) 6 , 
X (9) = -2.9 x 1(T 85 (m/V) f 



(18) 
(19) 
(20) 
(21) 



entering Eq. ^ which satisfactorily reproduces the 
TDSE data for r = 20 fs. However, changing the 
pulse duration destroys this agreement [see results for 
HOKE_FIT, green dash-dotted line in Fig. |] (a), (b)]. 
In contrast, comparison to the Standard models with ei- 
ther ionization rate (red line = MPI rate, blue = tunnel 
rate) shows reasonable agreement for the physically ex- 
tracted model parameters [see Eqs. (12) and (15)-(17)] 
for all pulse durations. 



An analysis of the physical mechanism responsible for 
the observed nonlinear refractive index saturation in the 
TDSE simulations is detailed in Fig. [2](c). Here, the re- 
sulting contributions to the nonlinear polarization from 
the bound ('BB' green line), continuum ('CC blue line) 
and bound-continuum ( ! BC red line) electronic wave 
function [see Eq. (|3|] are shown for r = 20 fs. The 
bound (BB) and bound-continuum (BC) contributions 
to P/vx(wrj) are both positively valued (focusing action) 
up to 45 TW/cm 2 , well beyond the threshold intensity 
where the overall nonlinear polarization (solid black line) 
changes sign. In contrast, the continuum contribution 
(CC) is close to zero for small intensities / < 20 TW/cm 2 , 
because no ionization occurs, and yields a negative (de- 
focusing) contribution when free charge generation kicks 
in. The value of this negative contribution (CC) yields 
significantly lower intensity thresholds than both the pos- 
itive contributions of BB and BC. Thus, the observed 
sign flip in the nonlinear polarization, i.e., the saturation 
of the nonlinear induced refractive index, can be clearly 
attributed to the contribution of free electrons according 
to Eq. ([3]). In this context it is worth mentioning that, 
strictly speaking, only the overall polarization is gauge 



G 



invariant. The projected contributions (BB, BC, CC) 
may differ for a different gauge in the TDSE equation. 
However, in particular the CC contribution preserves its 
qualitative focusing properties using different gauges, as 
shown in |13j . 

To quantify the impact of the free electrons, ioniza- 
tion is considered in Fig. [2] (d). The free electron 
density obtained from TDSE simulations via p(t) — 
p a t J dx\^ c (x:t)\ 2 [shown for I = 30 TW/cm 2 and r = 
20 fs; solid black line] develops strong oscillations at times 
where the pulse interacts with the atom. This is due to 
the fact that our definition using Eq. ([2| for the free elec- 
tron density is only applicable for the field free situation 
after the driving pulse has passed. In the presence of the 
laser field the eigenstates of H do no longer represent the 
bound states of the laser dressed atom. Therefore, spuri- 
ous oscillations occur due to contributions from electrons 
in high lying bound (Rydberg/Kramers-Henneberger) 
states close to the bulge of the bended atomic poten- 
tial and from truly free electrons. However, in the Stan- 
dard models we apply ionization rates that describe the 
truly ionized fraction after the laser pulse only, and the 
corresponding charge contributions do not feature these 
oscillations [red line = MPI rate, blue line = tunnel rate 
in Fig. [2] (d)]. Importantly, the truly ionized fraction 
according to the Standard models [dashed blue line for 
tunnel rate in [2] (c)] reproduces the CC contributions 
from TDSE simulations for all intensities and pulse du- 
rations. Minor deviations for / > 30 TW/cm 2 are due 
to previously discussed resonant occupation of high ly- 
ing Rydberg states [2U], whose description is beyond the 
scope of simple MPI/tunnel ionization models. 



B. Behavior in Temporal Domain 

Let us continue our comparison of TDSE results for 
Pnl with approximate models in time domain. In Fig. 
[3] (a) the time evolution of Pnl from TDSE simulations 
and the approximate HOKE_FIT and Standard models 
are shown for r = 20 fs and I = 20 TW/cm 2 . For 
that intensity, charge contributions are small and all 
three models reproduce accurately the response of the 
bound electrons, instantaneously following the electric 
field. For / = 35 TW/cm 2 in (b) charge contributions 
become important, resulting in a delayed response of op- 
posite sign kicking in at t ~ 5 fs, according to TDSE 
simulations (black solid line). The HOKEJTT model 
(dash-dotted green line), although yielding a compara- 
ble value of P/vl(^o), cannot describe that delayed con- 
tribution due to its instantaneous nature and produces 
much higher polarization amplitudes than observed in 
the TDSE simulations. In contrast, both Standard mod- 
els in (c), (d) qualitatively capture this delayed response 
and they develop only limited differences in amplitude. A 
similar situation is encountered in (e) for r = 100 fs and 
I = 24 TW/cm 2 , where also the limits of the Standard 
models are illustrated more clearly: The qualitative tem- 



poral character is described correctly over long evolution 
times (oscillation in phase with TDSE results). However, 
due to a slightly overestimated charge contribution, the 
Standard models yield larger amplitudes for the delayed 
response coming along with an inaccurate onset of the 
latter. 



IV. COMPARING TO 3D TDSE 
CALCULATIONS 

Finally, let us confirm the results obtained from the 
ID models by performing the same analysis in 3D. Here, 
we numerically solve 

it&tW, *) = (Ho + H mt ) t) (22) 
with H a = -V 2 /(2m e ) + V(r) including V(r) = 



-g 2 /47reor with 



\r\. The interaction with the 



linearly polarized external laser field in length gauge is 
computed in dipole approximation H int = —q e E(t)x |37j . 
We followed the same procedure to extract the param- 
eters of the approximate models as in the ID case (see 
Sec. pCl). 



From low intensity calculations we obtain a linear sus- 
ceptibility of x — 2.3 x 10~ 4 . The nonlinear suscep- 
tibilities x t X > ' ' ' evaluated at harmonic frequencies 
[see Eq. ( 10 )] are presented in Figs. E] (a)-(c) for pulse 
durations of r — 20 fs (solid black line) and r = 100 fs 
(solid magenta line). The deviations for different r are 
similar to those observed in the ID case, indicating simi- 
lar nonlinear dispersion. Because we again meet difficul- 
ties in extracting a consistent set of X s from Eqs. (Ill, 



we omit a repeated comparison to the HOKEJ35 model. 
Thus we resort to comparisons involving the HOKE_FIT 
model, defined through 



x ( 3 ) = 0.65 x 10~ 25 (m/V) 2 , 



,(5) 

X (7) = -1.0 x 10"°° (m/V) c 
X (9) = -1.0 x 10~ 87 (m/V) 6 



X'"' = 0.13 x 10~ 45 (m/V) 4 , 

-66 , 



and the Standard models employing 

x ( 3 ) = 0.65 x 10~ 25 (m/V) 2 

together with (K = 6) 

cr 6 = 1.13 x 10" 30 (cm 2 /TW) 6 /s 

for the MPI rate or 

a = 2.0 x 10 26 V/ms, 
/3 = 1.9 x 10 n V/m 



(23) 
(24) 
(25) 
(26) 

(27) 

(28) 



(29) 
(30) 



for the tunnel ionization rate. Naturally, the values of 
the model parameters in 3D differ from the ones obtained 
from ID TDSE data. Nevertheless, our 3D calculations 
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FIG. 3: ID nonlinear polarization in time domain, (a) Results from TDSE and approximate models for r = 20 fs at / = 
20 TW/cm 2 , where free electron contributions are negligible and the response is governed by the instantaneous contribution from 
bound electrons, (b)-(e) Temporal shapes, where delayed charge contributions become important: (b) TDSE and HOKE_FIT, 
(c) TDSE and Standard model with MPI rate and (d) TDSE and Standard model with tunnel rate for r = 20 fs and 
I = 35 TW/cm 2 ; (e) TDSE and MPI/tunnel Standard models for r = 100 fs and J = 24 TW/cm 2 . All temporal field shapes 
were filtered in Fourier domain using a hyperbolic tangent mask around 0.8ojo < CJ < 1.2wo- Color coding is the same as in 
Fig- i 



yield the same order of magnitude for the x^' s as ob- 
served in ID. It is interesting that x'' 5 '* an d X^ 7 "* exhibit 
opposite signs compared to the ID case. 

Let us continue with the comparison of 3D TDSE data 
with HOKE_FIT and Standard model results. In Fig. 
[5] (a) for t = 20 fs and (b) for r = 100 fs the nonlin- 
ear polarization Pnl{^o) inferred from 3D simulations 
(solid black line) is compared to the HOKE_FIT model 
(dash-dotted green line) and to the Standard models with 
a MPI rate (solid red line) or a tunnel ionization rate 
(solid blue line). As in the ID case, the HOKE_FIT 
model exhibits strong deviations for different pulse du- 
rations r. In contrast, the Standard models with physi- 
cally extracted parameters reasonably approximate sim- 
ulation results for both pulse durations r. However, for 
intensities I > 50 TW/cm 2 for the r = 20 fs case and 
I > 30 TW/cm 2 for r = 100 fs the overestimation of 
charge contributions seems more pronounced than in ID. 
In Fig. [5] (c) and (d) this is evident when comparing the 
temporal evolution of -Pjvx(i) from 3D TDSE simulations 
(solid black line) to Standard models with (c) MPI (red 
solid line) and (d) tunnel ionization rate (blue solid line) . 
Here, we show data for r = 20 fs and / = 50 TW/cm 2 , 
which ensures a considerable charge contribution. Both 



Standard models qualitatively reproduce simulation re- 
sults by capturing the instantaneous and delayed prop- 
erty. However, deviations in the amplitude of the delayed 
contribution are larger than in similar ID situations [Fig. 
|] (c), (d)]. The HOKE_FIT model (dash-dotted green 
line) in Fig. [Hj(e) for r = 100 fs and 1 = 40 TW/cm 2 not 
even qualitatively captures the delayed temporal prop- 
erty as it gives out of phase oscillations with an opposite 
sign for the field amplitude for t > 10 fs. 



V. CONCLUSION 

In conclusion, numerical solution of the TDSE for 
atomic hydrogen in ID and 3D reveals the saturation 
and change of sign of the nonlinear polarization at central 
laser frequency for rising intensity of the driving pulse. 
As in [13] . we identify a time retarded contribution re- 
sponsible for the transition from the focusing to defo- 
cusing regime. This time retarded contribution cannot 
be modeled by means of an instantaneous HOKE model, 
as we neither can physically extract consistent HOKE 
terms nor (even qualitatively) reproduce temporal field 
structures from simulations with a phenomenologically 
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FIG. 4: Model parameters extracted from 3D TDSE simulations, (a) x (3) , (b) X (5) , (c) X (?) for r = 20 is (black), r = 100 fs 
(magenta) as function of the peak intensity I. (d) The MPI cross section o?, extracted from 3D TDSE simulations [same color 
coding as in (a)]. 



fitted HOKE.FIT model. On the other hand, applying 
Standard models with either a MPI or tunnel ionization 
rate, clearly attributes the delayed response to the gen- 
eration of free electrons. With both standard models 
we achieve reasonable agreement with TDSE simulations 
for the nonlinear polarization at central laser frequency 
as well as for the temporal field structures. We propose a 
unique set of parameters, valid for all applied pulse dura- 
tions and intensities, to model the nonlinear polarization. 
This simple approach opens the possibility to efficiently 
simulate laser pulse propagation with reasonable accu- 



racy. However, due to their inherently approximate na- 
ture, possibly important details observed in TDSE simu- 
lations, e.g. the non-monotonous behavior of the ioniza- 
tion probability due to channel closure and resonant ex- 
citation, are not included in the Standard models. Thus, 
in order to accurately describe laser pulse propagation, a 
coupling of TDSE calculations to propagation equations 
as in |38| seems mandatory. 

We gratefully thank M. Ivanov, M. Kolesik, A.M. 
Popov and E. Cormier for fruitful and inspiring discus- 
sions on the bound electrons in the continuum. 
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